#include <GL/gl.h>

#include "Flob.h"

using std::swap;

Flob::Flob(float n, float d): sft_a(&Flob::t_threadFunc, &tf_a),
		sft_b(&Flob::t_threadFunc, &tf_b),
		sft_c(&Flob::t_threadFunc, &tf_c){
	inited=false;
	Nvalue=n;
	diff=d;
	velx=vely=velx0=vely0=accx=accy=dens=dens0=diverg=NULL;
	size=0;
	tf_c.myFlob=tf_b.myFlob=tf_a.myFlob=this;
	tf_c.funcToRun=tf_b.funcToRun=tf_a.funcToRun=TF_DATA::TF_NONE;
	tf_c.dt=tf_b.dt=tf_a.dt=0;
}
Flob::~Flob(){
	delete [] velx;
	delete [] vely;
	delete [] velx0;
	delete [] vely0;
	delete [] accx;
	delete [] accy;
	delete [] dens;
	delete [] dens0;
	delete [] diverg;
}
vec2i Flob::resize(vec2i sxy){
	sxy/=FLOB_DIVISOR;
	if(sxy.x<4) sxy.x=4;
	if(sxy.y<4) sxy.y=4;
	if(!inited){
		//easy
		sizexy=sxy;
		size=sizexy.x*sizexy.y;
		velx	= new float[size];
		vely	= new float[size];
		velx0	= new float[size];
		vely0	= new float[size];
		accx	= new float[size];
		accy	= new float[size];
		dens	= new float[size];
		dens0	= new float[size];
		diverg	= new float[size];

		memset(velx, 0, size);
		memset(vely, 0, size);
		memset(velx0, 0, size);
		memset(vely0, 0, size);
		memset(accx, 0, size);
		memset(accy, 0, size);
		memset(dens, 0, size);
		memset(dens0, 0, size);
		memset(diverg, 0, size);

		inited=true;
	}
	else{
		//clean up unneeded stuff
		delete [] diverg;
		delete [] velx;
		delete [] vely;
		delete [] accx;
		delete [] accy;
		delete [] dens;
		//set up
		size_t nsize=sxy.x*sxy.y;
		int minx = std::min(sxy.x, sizexy.x);
		int miny = std::min(sxy.y, sizexy.y);
		//preserve old
		velx = new float[nsize];
		vely = new float[nsize];
		dens = new float[nsize];
		accx = new float[nsize];
		accy = new float[nsize];
		diverg = new float[nsize];
		for(size_t i=0; i<nsize; ++i){
			velx[i] = vely[i] = dens[i] = accx[i] = accy[i] = diverg[i] = 0.f;
		}
		for(int y=0; y<miny; ++y)
			for(int x=0; x<minx; ++x){
				velx[y*sxy.x+x]=velx0[IX(x,y)];
				vely[y*sxy.x+x]=vely0[IX(x,y)];
				dens[y*sxy.x+x]=dens0[IX(x,y)];
			}
		//resize old
		delete [] velx0;
		delete [] vely0;
		delete [] dens0;
		velx0 = new float[nsize];
		vely0 = new float[nsize];
		dens0 = new float[nsize];
		memcpy(velx0, velx, nsize*sizeof(float));
		memcpy(vely0, vely, nsize*sizeof(float));
		memcpy(dens0, dens, nsize*sizeof(float));
		//finalize sizes
		size=nsize;
		sizexy=sxy;

	}
	dTrack.resize(size);
	vxTrack.resize(size);
	vyTrack.resize(size);
	return sizexy*FLOB_DIVISOR;
}
vec2i Flob::getSize() const{
	return sizexy;
}
////HEADACHE RELIEVERS////
float Flob::getVelx(vec2f pos) const{
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-2;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.y=sizexy.y-2;
	vec2i posi=pos;
	pos-=posi;
	size_t index=IX(posi);
	return (1-pos.y)*((1-pos.x)*(velx0[index])+(pos.x)*(velx0[index+1])) +
			(pos.y)*((1-pos.x)*(velx0[index+sizexy.x])+(pos.x)*(velx0[index+sizexy.x+1]));
}
float Flob::getVely(vec2f pos) const{
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-2;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.y=sizexy.y-2;
	vec2i posi=pos;
	pos-=posi;
	size_t index=IX(posi);
	return (1-pos.y)*((1-pos.x)*(vely0[index])+(pos.x)*(vely0[index+1])) +
			(pos.y)*((1-pos.x)*(vely0[index+sizexy.x])+(pos.x)*(vely0[index+sizexy.x+1]));
}
vec2f Flob::getVel(vec2f pos) const{
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-2;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.y=sizexy.y-2;
	vec2i posi=pos;
	pos-=posi;
	size_t index=IX(posi);
	return (1-pos.y)*((1-pos.x)*(vec2f(velx0[index], vely0[index]))+(pos.x)*(vec2f(velx0[index+1], vely0[index+1]))) +
			(pos.y)*((1-pos.x)*(vec2f(velx0[index+sizexy.x], vely0[index+sizexy.x]))+(pos.x)*(vec2f(velx0[index+sizexy.x+1], vely0[index+sizexy.x+1])));
}
float Flob::getDens(vec2f pos) const{
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-2;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.y=sizexy.y-2;
	vec2i posi=pos;
	pos-=posi;
	size_t index=IX(posi);
	return (1-pos.y)*((1-pos.x)*(dens0[index])+(pos.x)*(dens0[index+1])) +
			(pos.y)*((1-pos.x)*(dens0[index+sizexy.x])+(pos.x)*(dens0[index+sizexy.x+1]));
}
float Flob::getVal(vec2f pos, float* d) const{
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-2;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.y=sizexy.y-2;
	vec2i posi=pos;
	pos-=posi;
	size_t index=IX(posi);
	return (1-pos.y)*((1-pos.x)*(d[index])+(pos.x)*(d[index+1])) +
			(pos.y)*((1-pos.x)*(d[index+sizexy.x])+(pos.x)*(d[index+sizexy.x+1]));
}

void Flob::setDens(vec2f pos, float sDens){
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-2;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.y=sizexy.y-2;
	vec2i posi=pos;
	pos-=posi;
	size_t index=IX(posi);
	dens0[index]=(1-pos.y)*(1-pos.x)*sDens;
	dens0[index+1]=(1-pos.y)*(pos.x)*sDens;
	dens0[index+sizexy.x]=(pos.y)*(1-pos.x)*sDens;
	dens0[index+sizexy.x+1]=(pos.y)*(pos.x)*sDens;
}
void Flob::addDens(vec2f pos, float aDens){
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-2;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.y=sizexy.y-2;
	vec2i posi=pos;
	pos-=posi;
	size_t index=IX(posi);
	dens0[index]+=(1-pos.y)*(1-pos.x)*aDens;
	dens0[index+1]+=(1-pos.y)*(pos.x)*aDens;
	dens0[index+sizexy.x]+=(pos.y)*(1-pos.x)*aDens;
	dens0[index+sizexy.x+1]+=(pos.y)*(pos.x)*aDens;
}

void Flob::setDens(vec2i pos, float sDens){
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-1;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.x=sizexy.y-1;
	dens0[IX(pos)]=sDens;
}
void Flob::setDens(size_t index, float sDens){
	dens0[index]=sDens;
}
void Flob::addDens(vec2i pos, float aDens){
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-1;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.x=sizexy.y-1;
	dens0[IX(pos)]+=aDens;
}
void Flob::addDens(size_t index, float aDens){
	dens0[index]+=aDens;
}

void Flob::addAcc(vec2f pos, vec2f force){
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-2;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.y=sizexy.y-2;
	vec2i posi=pos;
	pos-=posi;
	size_t index=IX(posi);
	addAcc(index, force*(1-pos.y)*(1-pos.x));
	addAcc(index+1, force*(1-pos.y)*(pos.x));
	addAcc(index+sizexy.x, force*(pos.y)*(1-pos.x));
	addAcc(index+sizexy.x+1, force*(pos.y)*(pos.x));
}
void Flob::addAcc(vec2i pos, vec2f force){
	if(pos.x<0) pos.x=0;
	if(pos.x>sizexy.x-1) pos.x=sizexy.x-1;
	if(pos.y<0) pos.y=0;
	if(pos.y>sizexy.y-1) pos.x=sizexy.y-1;
	size_t index=IX(pos);
	accx[index]+=force.x;
	accy[index]+=force.y;
}
void Flob::addAcc(size_t index, vec2f force){
	accx[index]+=force.x;
	accy[index]+=force.y;
}

////HEADACHE GENERATORS////
void Flob::setBnd(int b, float* x, float def){
	if(b==0){
		for(int i=0; i<sizexy.x; ++i){
			x[i]=def;
			x[(sizexy.y-1)*sizexy.x+i]=def;
		}
		for(int i=0; i<sizexy.y; ++i){
			x[i*sizexy.x]=def;
			x[i*sizexy.x+sizexy.x-1]=def;
		}
	}
	//return;
	if(b==1){
		for(int i=0; i<sizexy.y; ++i){
			x[i*sizexy.x]=FLOB_BORDER_BOUNCE*x[i*sizexy.x+1];
			x[i*sizexy.x+sizexy.x-1]=FLOB_BORDER_BOUNCE*x[i*sizexy.x+sizexy.x-2];
		}
		for(int i=0; i<sizexy.x; ++i){
			x[i]=def;//x[i+sizexy.x];
			x[(sizexy.y-1)*sizexy.x+i]=def;//x[(sizexy.y-2)*sizexy.x+i];
		}
	}
	else if(b==2){
		for(int i=0; i<sizexy.x; ++i){
			x[i]=FLOB_BORDER_BOUNCE*x[i+sizexy.x];
			x[(sizexy.y-1)*sizexy.x+i]=FLOB_BORDER_BOUNCE*x[(sizexy.y-2)*sizexy.x+i];
		}
		for(int i=0; i<sizexy.y; ++i){
			x[i*sizexy.x]=def;//x[i*sizexy.x+1];
			x[i*sizexy.x+sizexy.x-1]=def;//x[i*sizexy.x+sizexy.x-2];
		}
	}
	else if(b==3){
		for(int i=0; i<sizexy.y; ++i){
			x[i*sizexy.x]=FLOB_BORDER_BOUNCE*x[i*sizexy.x+1];
			x[i*sizexy.x+sizexy.x-1]=FLOB_BORDER_BOUNCE*x[i*sizexy.x+sizexy.x-2];
		}
		for(int i=0; i<sizexy.x; ++i){
			x[i]=FLOB_BORDER_BOUNCE*x[i+sizexy.x];
			x[(sizexy.y-1)*sizexy.x+i]=FLOB_BORDER_BOUNCE*x[(sizexy.y-2)*sizexy.x+i];
		}
	}
	else if(b==4){
		for(int i=0; i<sizexy.y; ++i){
			x[i*sizexy.x]=x[i*sizexy.x+1];
			x[i*sizexy.x+sizexy.x-1]=x[i*sizexy.x+sizexy.x-2];
		}
		for(int i=0; i<sizexy.x; ++i){
			x[i]=x[i+sizexy.x];
			x[(sizexy.y-1)*sizexy.x+i]=x[(sizexy.y-2)*sizexy.x+i];
		}
	}
	//x[0]=x[size-1]=x[sizexy.x-1]=x[size-sizexy.x-1]=def;
}
void Flob::diffuseVal(int b, int iter, float rate, float* d, float* d0, float def){
	size_t index;
	float recip=1.0f/(1+4*rate);
	for(int k=0; k<iter; ++k){
		for(int y=1; y<sizexy.y-1; ++y)
			for(int x=1; x<sizexy.x-1; ++x){
				index=IX(x,y);
				d[index] = (d[index] + rate*(
						d0[index+1] +
						d0[index-1] +
						d0[index+sizexy.x] +
						d0[index-sizexy.x]))*recip;
			}
		setBnd(b, d, def);

	}
}
void Flob::advectVal(int b, float rate, float* d, float* d0, float def){
	size_t index;
	for(int y=1; y<sizexy.y-1; ++y)
		for(int x=1; x<sizexy.x-1; ++x){
			index=IX(x,y);
			d[index]=getVal(vec2f(x-velx0[index]*rate, y-vely0[index]*rate), d0);
		}
	setBnd(b, d, def);
}
void Flob::project(int iter, float dVel, float dDens, float dt){
	size_t index;
	float mult=0.5f/Nvalue;
	for(int y=1; y<sizexy.y-1; ++y)
		for(int x=1; x<sizexy.x-1; ++x){
			index=IX(x,y);
			diverg[index] = (velx0[index-1]-velx0[index+1] +
					vely0[index-sizexy.x]-vely0[index+sizexy.x])*mult;
			dens[index] = (diverg[index]+dens0[index-1]+dens0[index+1]+
					dens0[index-sizexy.x]+dens0[index+sizexy.x])/4.0f;
		}
	setBnd(0, dens, dDens);
	mult=0.5f*Nvalue;
	for(int y=1; y<sizexy.y-1; ++y)
		for(int x=1; x<sizexy.x-1; ++x){
			index=IX(x,y);
			velx[index]+=mult*(dens[index-1]-dens[index+1]);
			vely[index]+=mult*(dens[index-sizexy.x]-dens[index+sizexy.x]);
			velx[index]*=0.99f;
			vely[index]*=0.99f;
		}
	setBnd(1, velx, dVel);
	setBnd(2, vely, dVel);

}
void Flob::projectThread(int iter, float dDens, float dt){
	size_t index;
	float mult=0.5f/Nvalue;
	for(int y=1; y<sizexy.y-1; ++y)
		for(int x=1; x<sizexy.x-1; ++x){
			index=IX(x,y);
			diverg[index] = (velx0[index-1]-velx0[index+1] +
					vely0[index-sizexy.x]-vely0[index+sizexy.x])*mult;
			dens[index] = (diverg[index]+dens0[index-1]+dens0[index+1]+
					dens0[index-sizexy.x]+dens0[index+sizexy.x])/4.0f;
		}
	setBnd(0, dens, dDens);
	//threaded call to velx and vely projects
	tf_a.dt=tf_b.dt=dt;
	tf_a.funcToRun=TF_DATA::TF_PROJ2;
	tf_b.funcToRun=TF_DATA::TF_PROJ3;
	sft_a.Launch();
	sft_b.Launch();
	sft_a.Wait();
	sft_b.Wait();

}
void Flob::projectVely(){
	size_t index;
	float mult=0.5f*Nvalue;
	for(int y=1; y<sizexy.y-1; ++y){
		for(int x=1; x<sizexy.x-1; ++x){
			index=IX(x,y);
			vely[index]+=mult*(dens[index-sizexy.x]-dens[index+sizexy.x]);
			vely[index]*=0.99f;
		}
	}
	setBnd(2, vely, 0);
}
void Flob::projectVelx(){
	size_t index;
	float mult=0.5f*Nvalue;
	for(int y=1; y<sizexy.y-1; ++y){
		for(int x=1; x<sizexy.x-1; ++x){
			index=IX(x,y);
			velx[index]+=mult*(dens[index-1]-dens[index+1]);
			velx[index]*=0.99f;
		}
	}
	setBnd(1, velx, 0);
}

void Flob::simDensStep(float dt){
	float rate = dt*diff*Nvalue*Nvalue;
	diffuseVal(4, 1, rate, dens, dens0, 1.0f);
	swap(dens, dens0);
	swap(velx, velx0);
	swap(vely, vely0);
	diffuseVal(4, 1, rate, dens, dens0, 1.0f);
	swap(dens, dens0);
	swap(velx, velx0);
	swap(vely, vely0);
	advectVal(0, dt/FLOB_DIVISOR, dens, dens0, 1.0f);
	swap(dens, dens0);
}
void Flob::simVelStep(float dt){
	float rate = dt*diff*Nvalue*Nvalue;
	diffuseVal(1, 1, rate, velx, velx0);
	diffuseVal(2, 1, rate, vely, vely0);
	swap(velx, velx0);
	swap(vely, vely0);
	swap(dens, dens0);
	diffuseVal(1, 1, rate, velx, velx0);
	diffuseVal(2, 1, rate, vely, vely0);
	swap(velx, velx0);
	swap(vely, vely0);
	swap(dens, dens0);
	project(1,0,1,dt);
	swap(velx, velx0);
	swap(vely, vely0);
	swap(dens, dens0);
	rate=dt/FLOB_DIVISOR;
	advectVal(0, rate, velx, velx0);
	advectVal(0, rate, vely, vely0);
	swap(velx, velx0);
	swap(vely, vely0);
	swap(dens, dens0);
	project(1,0,1,dt);
	swap(velx, velx0);
	swap(vely, vely0);
	swap(dens, dens0);
}
void Flob::simFlStep(float dt){
	simDiffAll(dt);
	//at this point everything has been diffused twice
	projectThread(1,1,dt);
	swap(dens, dens0);
	swap(velx, velx0);
	swap(vely, vely0);
	simAdvAll(dt);
	swap(dens, dens0);
	swap(velx, velx0);
	swap(vely, vely0);
	tf_a.funcToRun=TF_DATA::TF_DISP_SCAN1;
	tf_b.funcToRun=TF_DATA::TF_DISP_SCAN2;
	tf_c.funcToRun=TF_DATA::TF_DISP_SCAN3;
	sft_a.Launch();
	sft_b.Launch();
	sft_c.Launch();
	sft_a.Wait();
	sft_b.Wait();
	sft_c.Wait();
}
void Flob::simPhysStep(float dt){
	tf_a.dt=tf_b.dt=dt;
	tf_a.funcToRun=TF_DATA::TF_ACC_VELX;
	tf_b.funcToRun=TF_DATA::TF_ACC_VELY;
	sft_a.Launch();
	sft_b.Launch();
	sft_a.Wait();
	sft_b.Wait();
}

void Flob::simDiffAll(float dt){
	tf_a.dt=tf_b.dt=tf_c.dt=dt;
	tf_a.funcToRun=TF_DATA::TF_DIFF_VELX;
	tf_b.funcToRun=TF_DATA::TF_DIFF_VELY;
	tf_c.funcToRun=TF_DATA::TF_DIFF_DENS;
	sft_a.Launch();
	sft_b.Launch();
	sft_c.Launch();
	sft_a.Wait();
	sft_b.Wait();
	sft_c.Wait();
	sft_a.Launch();
	sft_b.Launch();
	sft_c.Launch();
	sft_a.Wait();
	sft_b.Wait();
	sft_c.Wait();
}
void Flob::simAdvAll(float dt){
	tf_a.dt=tf_b.dt=tf_c.dt=dt;
	tf_a.funcToRun=TF_DATA::TF_ADV_VELX;
	tf_b.funcToRun=TF_DATA::TF_ADV_VELY;
	tf_c.funcToRun=TF_DATA::TF_ADV_DENS;
	sft_a.Launch();
	sft_b.Launch();
	sft_c.Launch();
	sft_a.Wait();
	sft_b.Wait();
	sft_c.Wait();
}

void Flob::drawDens(byte color){
	dTrack.doCalculations(dens0, size);
	float col[3] = {0.0f, 0.0f, 0.0f};
	size_t index=0;
	glBegin(GL_POINTS);
	for(int y=0; y<sizexy.y; ++y)
		for(int x=0; x<sizexy.x; ++x){
			col[color]=dTrack.normalTransform(dens0[IX(x,y)]);
			glColor3fv(col);
			glVertex2i(x,y);
		}
	glEnd();
}
void Flob::drawVel(){
	vxTrack.doCalculations(velx, size);
	vyTrack.doCalculations(vely, size);
	float col[3] = {0.0f, 0.0f, 0.0f};
	size_t index=0;
	glBegin(GL_POINTS);
	for(int y=0; y<sizexy.y; ++y)
		for(int x=0; x<sizexy.x; ++x){
			col[0]=vxTrack.normalTransform(velx[index]);
			col[2]=vyTrack.normalTransform(vely[index++]);
			glColor3fv(col);
			glVertex2i(x,y);
		}
	glEnd();
}
void Flob::drawSmart(){
	float col[3] = {0.0f, 0.0f, 0.0f};
	size_t index=0;
	glPushMatrix();
	glDisable(GL_TEXTURE_2D);
	glScalef(FLOB_DIVISOR, FLOB_DIVISOR, 1);
	glPointSize(ceil(FLOB_DIVISOR));
	glBegin(GL_POINTS);
	glColor4f(0.f,0.f,0.f,1.f);
	for(int y=0; y<sizexy.y; ++y)
		for(int x=0; x<sizexy.x; ++x){
			col[0]=vxTrack.normalTransform(fabs(velx[index]));
			col[1]=dTrack.normalTransform(dens[index]);
			col[2]=vyTrack.normalTransform(fabs(vely[index++]));
			glColor3fv(col);
			glVertex2i(x,y);
		}
	glEnd();
	glPointSize(1.0f);
	glPopMatrix();
}
void Flob::dbgPrintInfo(std::ostream& out){
	out << "Dens Avg:   " << dTrack.average() << '\n';
	out << "Dens Stdev: " << dTrack.deviation() << '\n';
	float velAvg = sqrt(vxTrack.average()*vxTrack.average()+vyTrack.average()*vyTrack.average());
	float velStdev = sqrt(vxTrack.deviation()*vxTrack.deviation()+vyTrack.deviation()*vyTrack.deviation());
	out << "Vel Avg:    " << velAvg << '\n';
	out << "Vel Stdev:  " << velStdev << '\n';
}

void Flob::getDensData(float& avg, float& stdev) const{
	avg=dTrack.average();
	stdev=dTrack.deviation();
}
void Flob::getVelxData(float& avg, float& stdev) const{
	avg=vxTrack.average();
	stdev=vxTrack.deviation();
}
void Flob::getVelyData(float& avg, float& stdev) const{
	avg=vyTrack.average();
	stdev=vyTrack.deviation();
}
void Flob::getVelData(float& avg, float& stdev) const{
	avg=sqrtf(vxTrack.average()*vxTrack.average()+vyTrack.average()*vyTrack.average());
	stdev=sqrtf(vxTrack.deviation()*vxTrack.deviation()+vyTrack.deviation()*vyTrack.deviation());
}

void Flob::t_threadFunc(void* dat){
	if(!dat) return;
	TF_DATA& tf = *static_cast<TF_DATA*>(dat);
	switch(tf.funcToRun){
	case TF_DATA::TF_NONE:
		break;
	case TF_DATA::TF_DIFF_DENS:{
		float rate = tf.dt*tf.myFlob->diff*tf.myFlob->Nvalue*tf.myFlob->Nvalue;
		tf.myFlob->diffuseVal(0, 1, rate, tf.myFlob->dens, tf.myFlob->dens0, 1.0f);
		swap(tf.myFlob->dens, tf.myFlob->dens0);
		break;
	}
	case TF_DATA::TF_DIFF_VELX:{
		float rate = tf.dt*tf.myFlob->diff*tf.myFlob->Nvalue*tf.myFlob->Nvalue;
		tf.myFlob->diffuseVal(4, 1, rate, tf.myFlob->velx, tf.myFlob->velx0, 1.0f);
		swap(tf.myFlob->velx, tf.myFlob->velx0);
		break;
	}
	case TF_DATA::TF_DIFF_VELY:{
		float rate = tf.dt*tf.myFlob->diff*tf.myFlob->Nvalue*tf.myFlob->Nvalue;
		tf.myFlob->diffuseVal(4, 1, rate, tf.myFlob->vely, tf.myFlob->vely0, 1.0f);
		swap(tf.myFlob->vely, tf.myFlob->vely0);
		break;
	}
	case TF_DATA::TF_ADV_DENS:{
		tf.myFlob->advectVal(0, tf.dt/FLOB_DIVISOR, tf.myFlob->dens, tf.myFlob->dens0, 1.0f);
		//swap(tf.myFlob->dens, tf.myFlob->dens0);
		break;
	}
	case TF_DATA::TF_ADV_VELX:{
		tf.myFlob->advectVal(4, tf.dt/FLOB_DIVISOR, tf.myFlob->velx, tf.myFlob->velx0, 1.0f);
		//swap(tf.myFlob->velx, tf.myFlob->velx0);
		break;
	}
	case TF_DATA::TF_ADV_VELY:{
		tf.myFlob->advectVal(4, tf.dt/FLOB_DIVISOR, tf.myFlob->vely, tf.myFlob->vely0, 1.0f);
		//swap(tf.myFlob->vely, tf.myFlob->vely0);
		break;
	}
	case TF_DATA::TF_ACC_VELX:{
		for(size_t i=0; i<tf.myFlob->size; ++i){
			tf.myFlob->velx[i]+=tf.myFlob->accx[i]*tf.dt;
			tf.myFlob->accx[i]=0;
		}
		break;
	}
	case TF_DATA::TF_ACC_VELY:{
		for(size_t i=0; i<tf.myFlob->size; ++i){
			tf.myFlob->vely[i]+=tf.myFlob->accy[i]*tf.dt;
			tf.myFlob->accy[i]=0;
		}
		break;
	}
	case TF_DATA::TF_PROJ1:
		break;
	case TF_DATA::TF_PROJ2:
		tf.myFlob->projectVelx();
		break;
	case TF_DATA::TF_PROJ3:
		tf.myFlob->projectVely();
		break;
	case TF_DATA::TF_DISP_SCAN1:
		tf.myFlob->vxTrack.doCalculations(tf.myFlob->velx, tf.myFlob->size, true);
		break;
	case TF_DATA::TF_DISP_SCAN2:
		tf.myFlob->vyTrack.doCalculations(tf.myFlob->vely, tf.myFlob->size, true);
		break;
	case TF_DATA::TF_DISP_SCAN3:
		tf.myFlob->dTrack.doCalculations(tf.myFlob->dens, tf.myFlob->size);
		break;
	case TF_DATA::TF_DISP_GEN:
		break;
	}
}

